Differential Inclusion Solver

Differential inclusions represent an important generalization of differential equations. The solution to a DI is a reachable set, instead of a single trajectory. These concepts are well known in control theory, but not yet exploited in the field of modeling and simulation. The main applications may be related to the uncertainty in dynamic systems. The solving procedure for differential inclusions is quite complicated compared to the numerical methods for differential equations. The algorithm is based on some concepts of the optimal con trol theory. Possible applications of differential inclusions in data uncertainty treatment are pointed out. It is shown that the differential inclusion is the most appropriate tool in dynamic uncertainty treatment. THE CHALLENGE There exists a s trange conviction among many simulationists who deal with continuous simulation that everything that happens in the real world can be described by differential equations. This approach to modeling is somewhat dangerous, because it forces the modellers to l ook for something (differential equation model) that might not exist at all. This may result in attempts to change the real world to fit our simulation tools, while the correct way should be completel opposite. The challenge is to look for new tools o to use those that have been known for a long time but not used or simply forgotten. Roughly speaking, a differential inclusion is similar to a differential equation, but its right-hand side is a set instead of a single value. While working with DIs we face another challenge. There are many theoretical results in this field, but little has been done in the area of numerical methods and practical applications. The solution to a DI needs hundreds of single system trajectory evaluations, which makes the whole task computationall expensive. Another challenge is the results representation (N-dimentional sets and the boundary surfaces). One of the main points of this paper is that ordinary (ODE) and partial differential equations as modeling tools are too primitive to describe the behavior of many real dynamical systems. The number of equations and internal feedbacks, nonlinearities and interactions are frequently used as a "measure" of the model complexity, while this should rather be related to the nature of the model and to the problem statement. APPLICATIONS IN DYNAMIC UNCERTAINTY The lack of reliable data in computer simulation is an important obstacle in many simulation projects. Models that are nice and valid from the academic point of view often result to be useless in practical applications, when the user cannot provide reliable data. In the past, a common way to treat this lack of exact data has been to suppose some model parameters or input variables to be random ones. This results in a stochastic model, where every realization of the system trajectory is different, and the problem is to determine the probability densit function in the system space for certain time sections, the variance, confidence intervals etc. Such stochastic anal ysis is interesting but not always possible. Some parameters of the model have "uncertain" values, and the model user has no idea about their probabilistic behavior. More likel we are given an interval that the uncertain parameter belongs to, instead of i ts probabilit distribution or sampled real data. Some external disturbances can fluctuate within certain intervals, and what we are asked to to give the interval for some output variables. The stochastic simulations with randomized variables do not give s uch intervals. Moreover, frequently the user wants to know a possible extremal value rather than a probability to reach them (recall the law of Murphy!). The uncertainty treatment has nothing, or very little, to do with "Monte Carlo" or stochastic simulation. The intervals we are looking for are not confidence intervals or any other statistics. There is no room here to mention a huge numbe of publication on uncertainty problems. See Bargiela (1998) for an example of uncertainty management in water distribution systems, for example. Another possible approach to uncertainty treatment are fuzzy sets. For more detail on the theory and applications of fuzzy sets consult Zadeh, Klir and Folger (see the References section). DIFFERENTIAL INCLUSION Let us consider a simple example of a second order system (1) dy/dt + a dy/dt + y = b This is a simple ODE model. Introducing notation x1= y , x2 = dy/dt we obtain (2) dx1/dt = x2 dx2/dt = b ax2 x1 In more general notation the state equation is (3) dx/dt = f(a,b,x) where x is a two-dimensional vector, t is the time and f is a vecto -valued function. Now suppose that the parameters a and b are uncertain and that the only information we have are the corresponding intervals where their values may belong, or a permissible (may be quite irregular and variable) set on the plain where the point (a,b) must belong. Note that we know nothing about a possible probability distribution of these parameters and we do not treat them as random variables. Thus, the equation (3) takes the following form. (4) dx/dt ∈ F(x,t) where F is a set. What we obtained is a differential inclusion (DI). The right-hand side of the equation (3) dete rmines the set F. However, this is merely one of several possible ways to represent the set. In this case it is parameterized b a and b. In general, the DI (4) should not necessarily be defined by the conditions (1), (2) and (3). The solution to a DI is the reachable set for the possible system trajectories that is exactly the solution to our uncertainty problem. In this very natural way we see that the uncertainty in dynamic system modeling leads to differential inclusions as a corresponding mathematical tool. This tool is known for about 70 years and that there is wide literature available on the DIs theory and applications. The first works have been published in 1931-32 by Marchaud and Zaremba. They used the terms "contingent" or "paratingent" eq uations. Later, in 1960-70, T. Wazewski and his collaborators published a series of works, referring to the DIs as orientor conditions and orientor fields. As always occurs with new theories, their works received sever criticism, mainly from some physicists who claimed that it is a stupid way of wasting time while dealing with so abstract an useless theories. Fortunately, the authors did not abandon the idea and developed the elemental theory of differential inclusions. In the decade 1930-40 such problems as the existence and properties of the solutions to the DIs have been resolved in the finite-dimensional space. After this, many works appear on DIs in more abstract, infinite-dimensional spaces. Within few years afte the first publications, the DIs resu lted to be the basic tool in the optimal control theory. Recall that optimal trajectories of a dynamic system are those that lay on the boundary of the system reachable set. In the works of Pontriagin, Lee and Markus, Bellman and many others, one of the f undamental problems is the determination of the properties of the reachable sets. Using the theory of Marchaud and Zaremba, T.Wazewski pointed out that in man optimal control problems the resulting control strategy is the so-called bang-bang control, generated by switching controllers. A more extended survey can be found in Raczynski, 1996. One of the best texts on DIs in finite-dimensional as well as abstract spaces is the book of Aubin and Cellina (1984). DIFFERENTIAL INCLUSION SOLVER A differential inclusion is a generalization of an ordinary differential equation. In fact, an ODE is a special case of a DI, where the right-hand F is a one-point set. One could expect that a solution algorithm for a DI might be obtained as some extension of known algorithms for the ODEs. Unfortunately, this is not the case. First of all, the solution to a DI is a set. Namely, it is a set in the time-state space, where the graphs of all possible trajectories of a DI are included. Finding the boundary of such set (named reachable set, or emission zone as in the works of Zaremba and Wazewski) is not an easy task. One of the properties of the reachable sets is the fact that if a trajectory reaches a point on the boundary of the RS at the final time, then its entire graph must belong to the RS. This fact is well known and used in the optimal control theory. Observe that any trajectory that reaches a point on the boundary of the RS is optimal in some sense. Such trajectories can be calculated using several methods, the main one being the Maximum Principle of Pontriagin . A similar method can be used to construct an algorithm for RS determination. If we can calculate a sufficient number of trajectories that scan the RS boundary, then we can see its shape. The trajectories should be uniformly distributed over the RS boundary. This can be done by some kind of random shooting over the RS boundary. Such shooting has nothing to do with a simple or primitive random shooting, when the trajectories are generated randomly inside the RS. My first attempts to develop a DI solver were presented on the IFAC Symposium on Optimization Methods, Varna, 1974. This was a random shooting method, but not a simple shooting. That algorithm generated trajectories inside the RS, but the control variable was being modified to obtain a nearl uniform distribution of points inside the RS at the end of the simulated time interval. The DI solve presented here is much more effective. In few words, the DI solver works as follows. The user provides the DI in the form of an equivalent control system. To do it he/she must parameterize the right-hand size (the set F) using a m-dimensional auxiliary variable u. The DI solver automatically generates the equations of the so called conjugate vecto p and integrates a set of trajectories, each of them belonging to the boundary of the RS. To achieve this, over each trajectory the Hamiltonian H(x,p,u,t) is maximized. This procedure is similar to that used in dynamic optimization. In the optimal control problem the main difficulty consists in the boundary conditions for the state and conjugate vectors. For the state vector we have the initial condition g iven, and for the conjugate vector only the final conditions (at the end of the trajectory) are known, given by the transversality conditions. This means that the optimal control algorithm must resolve the corresponding two-point-boundary value problem. In the case of a DI we are in a better situation. There is no objective function and no transversalit conditions. As the consequence, for the vecto p we can define the final as well as the initial conditions. Anyway we obtain a trajectory which graph belongs to the RS boundary. Defining initial conditions for p we integrate the trajectory only once, going forward. The only problem is how to generate these initial conditions in order to scan the RS boundary with nearly uniform density. The algorithm is quite simple: the initial conditions for p are generated randomly, due to a density function that is being automatically modified, covering with more density points that correspond to trajectories that fall into low density places at the RS boundary. Traject ries that are very close to each other are not stored (storing only one from each eventual cluster). As the result we obtain a set of trajectories covering the RS boundary that can be observed in graphical form and processed. EXAMPLE Consider a second order dynamic system where both the acceleration as well as the damping coefficient are uncertain. An example of the corresponding DI in parameterized form can be as follows. dx1/dt = x2 dx 2/dt = u1 u2x2 x1 where -1 < u1 < 1 and 0.5 < u2 < 1.5 . Figure 1 shows the 3D image of the RS in coordinates x1 (upwards), x2 (to the right) and t. On Figure 2 you can see a time-section of the RS for some fixed time. Observe the small cluster of trajectories at the origin of the coordinate system. These are trajectories (10000 in total) obtained by a simple random shooting, where both controls had been generated randomly within the limits defined above. The other pixels (the RS contour) are the end points (400 in total) of the trajectories generated by the DI solver. This demonstrates yhe uselessness is the simple shooting method in DI solving. CONCLUSIONS Differential inclusions can and should be used in continuous simulation, because there are problems that can hardly be solved using ODE or partial differential equations. The challenge in continuous simulation is to look for more appropriate modeling tools. Solving DIs is itself a challenge, taking into account the computational difficulties and the complexity of corresponding numerical algorithms. Also note that the simple random shooting, used sometimes in uncertainty problems can result in completely wrong assessments of the solution to a differential inclusion.